#!/bin/tcsh -ef
#############################################################################
#                                                                           #
# Title: Refine Lattice                                                     #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 30/04/2008	                                            #
# Last Modification: 30/04/2008	                                            #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 57
#
# MANUAL: This script refines the lattice, using the MRC program mmlatref.
#
# DISPLAY: imagename
# DISPLAY: nonmaskimagename
# DISPLAY: imagenumber
# DISPLAY: imagesidelength
# DISPLAY: lattice
# DISPLAY: secondlattice
# DISPLAY: backuplattice
# DISPLAY: comment
# DISPLAY: Refine_Lattice
# DISPLAY: det_tilt_lat
# DISPLAY: TLTAXIS
# DISPLAY: TLTANG
# DISPLAY: TLTAXA
# DISPLAY: TAXA
# DISPLAY: TANGL
# DISPLAY: realang
# DISPLAY: realcell
# DISPLAY: SYM
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set imagename = ""
set nonmaskimagename = ""
set imagenumber = ""
set imagesidelength = ""
set lattice = ""
set Refine_Lattice = ""
set realcell = ""
set TLTAXIS = ""
set TLTANG = ""
set secondlattice = ""
set magnification = ""
set stepdigitizer = ""
set tempkeep = ""
set revhk = ""
set realang = ""
set RESMIN = ""
set RESMAX = ""
set ALAT = "" 
set det_tilt_lat = ""
set delta = ""
set peakNum = ""
set testHand = ""
#
#$end_vars
#
set scriptname = 2dx_refineLattice
#
\rm -f LOGS/${scriptname}.results
#
set IS_2DX = yes
source ${proc_2dx}/initialize
#
echo "<<@evaluate>>"
#
if ( ${Refine_Lattice} == "0" ) then
  ${proc_2dx}/linblock "Not running."
  exit
endif
#
set halimsid = ${imagesidelength}
@ halimsid /= 2
echo halimsid = ${halimsid}
#
echo "<<@progress: 20>>"
#
if ( ! -e FFTIR/${imagename}_fft.mrc ) then
  ${proc_2dx}/linblock "ERROR: FFTIR/${imagename}_fft.mrc does not exist."
  ${proc_2dx}/protest "ERROR: First run Calculate FFT."
endif
#
#
#################################################################################
${proc_2dx}/linblock "2dx_mmlatref - to refine the lattice."
#################################################################################
#
\rm -f TMP456111.tmp
#
echo "Running:"
cat << eot
${bin_2dx}/2dx_mmlatref.exe
FFTIR/${imagename}_fft.mrc
${imagenumber} ${nonmaskimagename} 
Y                       ! use grid units?
Y                       ! generate grid from lattice?
N                       ! generate points from lattice?
2 0 0 30 30 15 15       ! IPIXEL,IOUT,NUMSPOT,NOH,NOK,NHOR,NVERT
20 ${halimsid} ${halimsid} ${halimsid}     ! RINNER,ROUTER,XORIG,YORIG
100 80.0 0               ! NCYCLES,FSHIFT,ILIST for lattice param refn
${lattice}                 ! AX,AY,BX,BY
eot
#
echo " "
echo "Now starting:"
echo " "
#
${bin_2dx}/2dx_mmlatref.exe << eot
FFTIR/${imagename}_fft.mrc
${imagenumber} ${nonmaskimagename} 
Y                       ! use grid units?
Y                       ! generate grid from lattice?
N                       ! generate points from lattice?
2 0 0 30 30 15 15       ! IPIXEL,IOUT,NUMSPOT,NOH,NOK,NHOR,NVERT
20 ${halimsid} ${halimsid} ${halimsid}     ! RINNER,ROUTER,XORIG,YORIG
100 80.0 0               ! NCYCLES,FSHIFT,ILIST for lattice param refn
${lattice}                 ! AX,AY,BX,BY
eot
#
if ( ! -e TMP456111.tmp ) then
  ${proc_2dx}/protest "ERROR in 2dx_mmlatref.exe occured."
endif
#
set newlat = `cat TMP456111.tmp`
set lattice = `echo ${newlat} | sed 's/ /,/g' `
#
echo "set lattice = ${lattice}" >> LOGS/${scriptname}.results
# echo "set Refine_Lattice = 0" >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "2dx_lencalc: To calculate some values for the tilt geometry determination."
#############################################################################
#
set lencalc_docfile = "2dx_lencalc-doc.tmp"
#
\rm -f ${lencalc_docfile}
#
${bin_2dx}/2dx_lencalc.exe << eot
${lattice}
${realcell}
${realang}
${imagesidelength}
${magnification}
${stepdigitizer}
${lencalc_docfile}
eot
#
cat ${lencalc_docfile}
#
set hand     = `head -n 2 ${lencalc_docfile} | tail -n 1`
set untilted = `head -n 4 ${lencalc_docfile} | tail -n 1`
set alpha    = `head -n 6 ${lencalc_docfile} | tail -n 1`
set tilted   = `head -n 8 ${lencalc_docfile} | tail -n 1`
#
\rm -f ${lencalc_docfile}
#
echo "<<@progress: 60>>"
#
#############################################################################
#                                                                           #
#  EMTILT                                                                   #
#                                                                           #
#                 ->  1   , 1   , 90.0       (or also : lenA,lenB,90.0)     #
#                 ->  lenH, lenK, ang(H,K)                                  #
#                                                                           #
#                 <-  ang(Tiltaxis to A*) = Omega                           #
#                 <-  Tiltangle (without sign)                              #
#                                                                           #
#  A* was in the EM in the sample-plane.                                    #
#  Omega is NOT the angle from the Tiltaxis to H, but this one without      #
#     perspective. Omega is in the negative plane.                          #
#     Omega is slightly smaller than the latter (TLTAXA).                   #
#                                                                           #
#############################################################################
#                                                                           #
#  Angle from TLT-axis to H is TLTAXA, TLT-axis to K is TLTAXB. (???)       #
#  Angle from Y-axis to H is TLTAXA, Y-axis to K is TLTAXB.                 #
#  Henderson: "TLTANG, TLTAXA, TLTAXB are on negative, nobody wants them."  #
#  This is not fully correct: TTMASK or TTBOX need them.                    #
#                                                                           #
#  TLTAXIS: Angle from X-axis to Tiltaxis                                   #
#  TLTANG is positive for strong underfocus at top of image.                #
#  TLTANG is positive for left overfocus, right underfocus when TLTAXIS= Y  #
#  TLTANG is positive for strong underfocus on left side, when TLTAXIS= Y   #
#       (Left screw positive counting)                                      #
#          TLTAXIS, TLTANG are for TTREFINE.                                #
#                                                                           #
#  TTREFINE has perfect output for ORIGTILTC :                              #
#          TAXA  , TAXB  , TANGL  are on sample, for ORIGTILTC.             #
#                                                                           #
#############################################################################
#
# Henderson in an email 6/2000:
#   y=0 means the x-axis which is horizontal.  It is indeed the bottom line of
# the image.  The first point on that line is (0,0), the second point (1,0)
# and so on [e.g. (x,0)].  Therefore all points on the line have y=0.  If the
# tiltaxis is horizontal, parallel to x and the defocus gets stronger as you
# go to higher values of y (i.e. further up the image), then TLTANGL is positive.
# This is a robust definition and TLTAXIS can take any number between -89.999
# and +89.999.  If it happens to be exactly 90.000 (note that this has never
# occurred in practice), then you have to know how the program reacts and that is
# where the extra x=0 description arises.  In practice, if your film orientation
# is such that the tiltaxis is mostly vertical but varies a bit with some of
# them on different sides of the vertical, then those with TLTAXIS of 85 degrees,
# for example, would have TLTANG positive, whereas those with TLTAXIS of 95 (i.e.
# -85) degrees would have TLTANGL negative.  So, the TLTANGL sign would change
# at 90.   This does not happen with the tiltaxis roughly horizontal, when there
# is no change in the definition as the TLTAXIS passes from -10 to +10.  I hope
#?| this is now clearer.
#
#############################################################################
#
#############################################################################
${proc_2dx}/linblock "emtilt: To calculate the tilt geometry from lattice distortions."
#############################################################################
#
\rm -f TMP_2dx_emtilt_1.out
#
set oldTLTAXIS = ${TLTAXIS}
set oldTLTANG  = ${TLTANG}
set emtilt_docfile = "2dx_emtilt-doc.tmp"
\rm -f ${emtilt_docfile}
#
echo untilted = ${untilted}
echo tilted = ${tilted}
echo oldTLTAXIS = ${oldTLTAXIS}
echo oldTLTANG = ${oldTLTANG}
echo hand = ${hand}
echo alpha = ${alpha}
#
${bin_2dx}/2dx_emtilt.exe << eot
${untilted}
${tilted}
${oldTLTAXIS}
${oldTLTANG}
${hand}
${alpha}
${emtilt_docfile}
eot
#
echo "2dx_emtilt.exe finished."
#
# read those values for the tilt geometry back into this script shell:
set LATTICE_TLTAXIS = `head -n 2  ${emtilt_docfile} | tail -n 1`
set LATTICE_TLTANG  = `head -n 4  ${emtilt_docfile} | tail -n 1`
set LATTICE_TLTAXA  = `head -n 6  ${emtilt_docfile} | tail -n 1`
set LATTICE_TAXA    = `head -n 8  ${emtilt_docfile} | tail -n 1`
set LATTICE_TANGL   = `head -n 10 ${emtilt_docfile} | tail -n 1`
\rm -f ${emtilt_docfile}
#
echo "<<@progress: 75>>"
#
set highangle = `echo ${LATTICE_TLTANG} | awk '{if($1>25 || $1<-25) {s=1} else {s=0}} END {print s}'`
#
if ( ${det_tilt_lat} == 'n' ) then
  set highangle = 0
endif
#
if ( ${highangle} == '1' ) then
  #############################################################################
  ${proc_2dx}/linblock "Saving tilt geometry from lattice, because of high tilt."
  #############################################################################
  set TLTAXIS = ${LATTICE_TLTAXIS}
  set TLTANG  = ${LATTICE_TLTANG}
  set TLTAXA  = ${LATTICE_TLTAXA}
  set TAXA    = ${LATTICE_TAXA}
  set TANGL   = ${LATTICE_TANGL}
  #
  # Make sure these values are going back into 2dx:
  echo "set TLTAXIS = ${TLTAXIS}" >> LOGS/${scriptname}.results
  echo "set TLTANG  = ${TLTANG}"  >> LOGS/${scriptname}.results
  echo "set TLTAXA  = ${TLTAXA}"  >> LOGS/${scriptname}.results
  echo "set TAXA    = ${TAXA}"    >> LOGS/${scriptname}.results
  echo "set TANGL   = ${TANGL}"   >> LOGS/${scriptname}.results
  echo "set DEFOCUS_ACTIVE = 2"   >> LOGS/${scriptname}.results
  #
else
  #############################################################################
  ${proc_2dx}/linblock "NOT saving tilt geometry from lattice, because of low tilt."
  #############################################################################
  echo "TLTAXIS = ${LATTICE_TLTAXIS}"
  echo "TLTANG  = ${LATTICE_TLTANG}" 
  echo "TLTAXA  = ${LATTICE_TLTAXA}" 
  echo "TAXA    = ${LATTICE_TAXA}"   
  echo "TANGL   = ${LATTICE_TANGL}"  
  echo "(These values were calculated, but will not be used.)"
endif
#
echo "set LATTICE_TLTAXIS = "\"${LATTICE_TLTAXIS}\" >> LOGS/${scriptname}.results
echo "set LATTICE_TLTANG  = "\"${LATTICE_TLTANG}\"  >> LOGS/${scriptname}.results
echo "set LATTICE_TLTAXA  = "\"${LATTICE_TLTAXA}\"  >> LOGS/${scriptname}.results
echo "set LATTICE_TAXA    = "\"${LATTICE_TAXA}\"    >> LOGS/${scriptname}.results
echo "set LATTICE_TANGL   = "\"${LATTICE_TANGL}\"   >> LOGS/${scriptname}.results
#
echo "# IMAGE-IMPORTANT: FFTIR/${imagename}_fft.mrc <FFT of Image>" >> LOGS/${scriptname}.results
echo "# IMAGE-IMPORTANT: FFTIR/${imagename}_red_fft.mrc <FFT of Downsampled Image>" >> LOGS/${scriptname}.results
#
echo " " >> History.dat
echo "::From Lattice:" >> History.dat
echo "::TLTAXIS = ${LATTICE_TLTAXIS}" >> History.dat
echo "::TLTANG  = ${LATTICE_TLTANG}" >> History.dat
echo "::TLTAXA  = ${LATTICE_TLTAXA}" >> History.dat
echo "::TAXA    = ${LATTICE_TAXA}" >> History.dat
echo "::TANGL   = ${LATTICE_TANGL}" >> History.dat
#
#############################################################################
${proc_2dx}/linblock "2dx_calcmag - to calculate the theoretical magnification"
#############################################################################
echo "<<@progress: 85>>"
echo "<<@evaluate>>"
#
set outputfile = 2dx_verifyParams.tmp
\rm -f ${outputfile}
#
${bin_2dx}/2dx_calcmag.exe << eot
${realcell}
${realang}
${TLTAXIS}
${TLTANG}
${lattice}
${imagesidelength}
${magnification}
${stepdigitizer}
${outputfile}
eot
#
if ( ! -e ${outputfile} ) then
  ${proc_2dx}/protest "ERROR in 2dx_calcmag.exe"
endif
#
set theormag = `cat ${outputfile} | head -n 1`
set RANGrec  = `cat ${outputfile} | head -n 2 | tail -n 1`
set RANGreal = `cat ${outputfile} | head -n 3 | tail -n 1`
#
\rm -f ${outputfile}
#
${proc_2dx}/linblock "Theoretical magnification is ${theormag}, given magnification is ${magnification}"
${proc_2dx}/linblock "Angle in reciprocal lattice is ${RANGrec}."
${proc_2dx}/linblock "Angle in real-space lattice is ${RANGreal}."
#
echo "set CALCULATEDMAG = ${theormag}" >> LOGS/${scriptname}.results
#
set docfile = 'peaks_xy_final.dat'
if ( -e ${docfile} ) then
  #############################################################################
  #                                                                           #
  ${proc_2dx}/linblock "2dx_laterror to find fitness of current lattice."
  #                                                                           #
  ############################################################################# 
  #
  ${bin_2dx}/2dx_laterror.exe "${lattice}" ${docfile}
  #
else
  #
  ${proc_2dx}/linblock "2dx_laterror not run, docfile with peaks not found."
  #
endif
#
echo "set LATTICE_done = y" >> LOGS/${scriptname}.results
echo "set SPOTS_done = n" >> LOGS/${scriptname}.results
echo "set UNBENDING_done = n" >> LOGS/${scriptname}.results
echo "set ML_done = n" >> LOGS/${scriptname}.results
echo "set CTF_done = n" >> LOGS/${scriptname}.results
echo "set MERGING_done = n" >> LOGS/${scriptname}.results
#
echo "<<@progress: 100>>"
#
##########################################################################
${proc_2dx}/linblock "${scriptname} - normal end."
##########################################################################
#
